//============================================================================
// Name        : csiena
// Author      : Josh Lospinoso
// Copyright   : MIT License
// Description : Austere, limited implementation of the R package RSiena in
//				 C++. Used for fast maximum likelihood and method of moments
//				 estimation of stochastic actor oriented models (SAOMs)
//				 without any frills or user experience niceties.
//============================================================================

#include "crimpEvaluation.h"

namespace csiena {

CrimpEvaluation::CrimpEvaluation(int actors) {
	this->actors = actors;
	#if GAMMA_LENGTH > 0
		this->gamma = new double[GAMMA_LENGTH];
		this->statistics = new double[2*GAMMA_LENGTH];  /* Needs to be revisited if there are more elaborate models! */
		this->score = new double[GAMMA_LENGTH];
		this->information = new double[GAMMA_LENGTH*GAMMA_LENGTH];
	#endif
	this->resetStatistics();
}

CrimpEvaluation::CrimpEvaluation(const CrimpEvaluation &donorEvaluation) {
	this->actors = donorEvaluation.actors;
	#if GAMMA_LENGTH > 0
		this->gamma = new double[GAMMA_LENGTH];
		this->statistics = new double[GAMMA_LENGTH];
		this->score = new double[GAMMA_LENGTH];
		this->information = new double[GAMMA_LENGTH];
	#endif
	this->logLikelihood = donorEvaluation.logLikelihood;
	for (int i=0; i<GAMMA_LENGTH; i++) {
		this->gamma[i] = donorEvaluation.gamma[i];
		this->score[i] = donorEvaluation.score[i];
		this->statistics[i] = donorEvaluation.statistics[i];
		for (int j=0; j<GAMMA_LENGTH; j++) {
			this->information[i*GAMMA_LENGTH+j] = donorEvaluation.information[i*GAMMA_LENGTH+j];
		}
	}

}

CrimpEvaluation::~CrimpEvaluation() {
#if GAMMA_LENGTH > 0
	delete(this->gamma);
	delete(this->statistics);
	delete(this->score);
	delete(this->information);
#endif
}

bool CrimpEvaluation::isCrimpEvaluation() {
	return true;
}

double CrimpEvaluation::getGamma(int gammaIndex) {
	if (gammaIndex < GAMMA_LENGTH) {
		return this->gamma[gammaIndex];
	} else {
		std::cerr << "CrimpEvaluation::getGamma: Index out of bounds" << gammaIndex << std::endl;
		return 0.0;
	}
}

void CrimpEvaluation::setGamma(int gammaIndex, double value) {
	if (gammaIndex < GAMMA_LENGTH) {
		if ( (value > MAX_GAMMA) | (value < MIN_GAMMA)) {
			std::cerr << "CrimpEvaluation::setGamma: Value out of bounds " << (value) << std::endl;
		} else {
			this->gamma[gammaIndex] = value;
		}
	} else {
		std::cerr << "CrimpEvaluation::setGamma: Index out of bounds" << gammaIndex << std::endl;
	}
}

double CrimpEvaluation::getLogLikelihood() {
	return this->logLikelihood;
}

double CrimpEvaluation::getScore(int gammaIndex) {
	if (gammaIndex < GAMMA_LENGTH) {
		return this->score[gammaIndex];
	} else {
		std::cerr << "CrimpEvaluation::getScore: Index out of bounds" << gammaIndex << std::endl;
		return 0.0;
	}
}

double CrimpEvaluation::getInformation(int gammaIndex1, int gammaIndex2) {
	if ((gammaIndex1 < GAMMA_LENGTH) & (gammaIndex2 < GAMMA_LENGTH)) {
		return this->information[gammaIndex1 * GAMMA_LENGTH + gammaIndex2];
	} else {
		std::cerr << "CrimpEvaluation::getInformation: Index out of bounds" << gammaIndex1 << "," << gammaIndex2 << std::endl;
		return 0.0;
	}
}

double CrimpEvaluation::getStatistics(int gammaIndex) {
	if (gammaIndex < GAMMA_LENGTH*2) {
		return this->statistics[gammaIndex];
	} else {
		std::cerr << "CrimpEvaluation::getStatistics: Index out of bounds" << gammaIndex << std::endl;
		return 0.0;
	}
}

void CrimpEvaluation::evaluate(Digraph* observedDigraph, Digraph* trueDigraph) {
	this->resetStatistics();
	if (GAMMA_LENGTH == 0) {
		// Classic model, so no informant accuracy. Just check that the two graphs are equal for debugging purposes.
		bool validCrimp=true;
		for (int i=0; i<this->actors;i++) {
			for (int j=0; j<this->actors; j++) {
				if (i != j) {
					if ( observedDigraph->get(i,j) != trueDigraph->get(i,j) ) {
						validCrimp=false;
						//break;
						std::cerr << std::endl<< i << " -> " << j << " dont match " << std::endl;
					}
				}
			}
		}
		if (!validCrimp) {
			std::cerr << "Error! the true digraph is not equal to the observed digraph and no informant accuracy model is specified." << std::endl;
		}
	} else if (GAMMA_LENGTH == 1) {
		// gamma(0) = error rate
		/* Statistic: | True: | Observed: |
		 *           0|      1|        1  |
		 *           1|      1|        0  |
		 *           0|      0|        0  |
		 *           1|      0|        1  |
		 */
		// Model 1, single error term. Count the differences as statistic 1
		for (int i=0; i<this->actors; i++) {
			for (int j=0; j<this->actors; j++) {
				if (i != j) {
					if ( observedDigraph->get(i,j) == trueDigraph->get(i,j) ) {
						this->statistics[0]++;
					} else {
						this->statistics[1]++;
					}
				}
			}
		}
		if ((this->getGamma(0) > 0) & this->getGamma(0) < 1) {
			this->logLikelihood = this->getStatistics(1) * log(this->getGamma(0)) + this->getStatistics(0) * log(1-this->getGamma(0));
			this->score[0] = this->getStatistics(1) / this->getGamma(0) - this->getStatistics(0) / (1-this->getGamma(0));
			this->information[0] = this->getStatistics(0) / ( (1-this->getGamma(0)) * (1-this->getGamma(0))) + this->getStatistics(1) / (this->getGamma(0)*this->getGamma(0));
		} else if ( (this->getGamma(0) == 0) | (this->getGamma(0) == 1) ) {
			//std::cerr << "Warning: with gamma value " << this->getGamma(0) << ", score and information are undefined." << std::endl;
			if (this->getGamma(0) == 0) {
				if (this->statistics[1] == 0) {
					this->logLikelihood = 0;
				} else {
					this->logLikelihood = -99999999;
				}
			} else {
				if (this->statistics[0] == 0) {
					this->logLikelihood = 0;
				} else {
					this->logLikelihood = -99999999;
				}
			}
		} else {
			std::cerr << "Error, gamma value out of range: " << this->getGamma(0) << std::endl;
		}
	} else if (GAMMA_LENGTH == 2) {
		// gamma(0) = false
		/* Statistic: | True: | Observed: | ToS notation:
		 *           0|      1|        1  | (1 - gamma(0) ) and g(0star)
		 *           1|      1|        0  | False negative, gamma(0) and g(0)
		 *           2|      0|        0  | (1 - gamma(1) ) and g(1star)
		 *           3|      0|        1  | False positive, gamma(1) and g(1)
		 */
		for (int i=0; i<this->actors; i++) {
			for (int j=0; j<this->actors; j++) {
				if (i != j) {
					if ( trueDigraph->get(i,j) ) {
						if (observedDigraph->get(i,j)) {
							this->statistics[0]++; // g(0star)
						} else {
							this->statistics[1]++; // g(0)
						}
					} else {
						if (!observedDigraph->get(i,j)) {
							this->statistics[2]++; // g(1star)
						} else {
							this->statistics[3]++; // g(1)
						}
					}
				}
			}
		}
		if ((this->getGamma(0) > 0) & this->getGamma(0) < 1) {
			this->logLikelihood =
					this->getStatistics(0) * log(1-this->getGamma(0)) +
					this->getStatistics(1) * log(this->getGamma(0)) +
					this->getStatistics(2) * log(1 - this->getGamma(1)) +
					this->getStatistics(3) * log(this->getGamma(1));
			this->score[0] = this->getStatistics(1) / this->getGamma(0) - this->getStatistics(0) / (1-this->getGamma(0));
			this->score[1] = this->getStatistics(3) / this->getGamma(1) - this->getStatistics(2) / (1-this->getGamma(1));
			this->setInformation(0,0, this->getStatistics(0) / ( (1-this->getGamma(0)) * (1-this->getGamma(0))) + this->getStatistics(1) / (this->getGamma(0)*this->getGamma(0)) );
			this->setInformation(1,1,  this->getStatistics(2) / ( (1-this->getGamma(1)) * (1-this->getGamma(1))) + this->getStatistics(3) / (this->getGamma(1)*this->getGamma(1)));
			this->setInformation(0,1, 0);
			this->setInformation(1,0, 0);
		} else if ( (this->getGamma(0) == 0) | (this->getGamma(0) == 1) ) {
			std::cerr << "Warning: with gamma value " << this->getGamma(0) << ", score and information are undefined." << std::endl;
			if (this->getGamma(0) == 0) {
				if (this->statistics[1] == 0) {
					this->logLikelihood = 0;
				} else {
					this->logLikelihood = -99999999;
				}
			} else {
				if (this->statistics[0] == 0) {
					this->logLikelihood = 0;
				} else {
					this->logLikelihood = -99999999;
				}
			}
		} else {
			std::cerr << "Error, gamma value out of range: " << this->getGamma(0) << std::endl;
		}
	}

}

void CrimpEvaluation::print() {
	std::cout << "Crimp evaluation:" << std::endl;
	if (GAMMA_LENGTH == 0) {
		std::cout << "  >> GAMMA=0, no informant accuracy model specified.\n" << std::endl;
	} else if (GAMMA_LENGTH ==1) {
		printf("  >> GAMMA=1: Single error rate specified.\n");
		printf("  >> Gamma = %+0.4f\n", this->getGamma(0));
		printf("  >> LogL. = %+0.4f\n", this->logLikelihood);
		printf("  >> Score = %+0.4f\n", this->getScore(0));
		printf("  >> Info. = %+0.4f\n", this->getInformation(0,0));
		printf("  >> Err's = %6u\n", (int)this->getStatistics(0));
	} else if (GAMMA_LENGTH ==2) {
		std::cout << "  >> GAMMA=2: Dual error rates specified.\n" << std::endl;
		printf("  >> Gamma(1->0) = %+0.4f\n", this->getGamma(0));
		printf("  >> Gamma(0->1) = %+0.4f\n", this->getGamma(1));
		printf("  >> LogL.       = %+0.4f\n", this->logLikelihood);
		printf("  >> Score(1->0) = %+0.4f\n", this->getScore(0));
		printf("  >> Score(0->1) = %+0.4f\n", this->getScore(1));
		printf("  >> Info.(1->0) = %+0.4f\n", this->getInformation(0,0));
		printf("  >> Info.(0->1) = %+0.4f\n", this->getInformation(1,1));
		printf("  >> #(1->1)     = %6u\n", (int)this->getStatistics(0));
		printf("  >> #(1->0)     = %6u\n", (int)this->getStatistics(1));
		printf("  >> #(0->0)     = %6u\n", (int)this->getStatistics(2));
		printf("  >> #(0->1)     = %6u\n", (int)this->getStatistics(3));
	} else {
		std::cerr << "Error, gamma value out of range: " << this->getGamma(0) << std::endl;
	}
}

void CrimpEvaluation::setStatistics(int gammaIndex, double value) {
	if (gammaIndex < GAMMA_LENGTH*2) {
		this->statistics[gammaIndex] = value;
	} else {
		std::cerr << "CrimpEvaluation::setStatistics: Index out of bounds" << gammaIndex << std::endl;
	}
}

void CrimpEvaluation::setScore(int gammaIndex, double value) {
	if (gammaIndex < GAMMA_LENGTH) {
		this->score[gammaIndex] = value;
	} else {
		std::cerr << "CrimpEvaluation::setScore: Index out of bounds" << gammaIndex << std::endl;
	}
}

void CrimpEvaluation::incrementGamma(int gammaIndex, double value) {
	if (gammaIndex < GAMMA_LENGTH) {
		this->setGamma(gammaIndex, this->getGamma(gammaIndex)+value);
	} else {
		std::cerr << "CrimpEvaluation::incrementGamma: Index out of bounds" << gammaIndex << std::endl;
	}
}

void CrimpEvaluation::setInformation(int gammaIndex1, int gammaIndex2, double value) {
	if ((gammaIndex1 < GAMMA_LENGTH) & (gammaIndex2 < GAMMA_LENGTH)) {
		this->information[gammaIndex1 * GAMMA_LENGTH + gammaIndex2] = value;
	} else {
		std::cerr << "CrimpEvaluation::setInformation: Index out of bounds" << gammaIndex1 <<"," << gammaIndex2 << std::endl;
	}

}

void CrimpEvaluation::resetStatistics() {
	for (int i=0; i<GAMMA_LENGTH; i++) {
		this->setScore(i, 0);

		for (int j=0; j<GAMMA_LENGTH; j++) {
			this->setInformation(i,j,0);
		}
	}
	for (int i=0; i<2*GAMMA_LENGTH; i++) {
		this->setStatistics(i,0);
	}
}

}
